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Abstract 

A self-consistent mean-field theory for bosons for T > is used to reconcile 
predictions of collapse with recent observations of Bose-Einstein condensa- 
tion of ^Li. Eigenfunctions of a (non-separable) Hamiltonian that includes 
the anisotropic external trap field and atom-atom interactions are obtained 
by an iteration process. A sum over the Bose distribution, and the "alternat- 
ing direction implicit" algorithm are used. Near Tc, the ensemble exhibits a 
localized condensate composed of atoms in the few lowest states. For lower 

T, numerical instability indicates collapse to a more dense phase. 
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Recent observations of Bose-Einstein condensation in tenuous gases of alkali atoms 
have brought to concrete realization a text-book paradigm and a prediction by Einstein 70 
years ago. These outstanding experimental achievements, coming after the rapid develop- 
ment of methods for cooling atoms by laser light and by evaporation in a trap, make it 
possible to test theories for weakly-interacting Bose gases. In contrast with superfluid he- 
lium, where the atoms interact strongly and thus depart from simple predictions, there is 
hope that BEC in dilute alkali atom gases can be usefully addressed by accessible theoretical 
methods. 

Perhaps the experimental result that presents the greatest challenge to understanding is 
the reported Bose-Einstein condensation of ^Li, which is known to have a scattering length, 
a = — 27.2ao @], indicating an attractive atom-atom interaction. It has been shown that in 
free space, BEC with a < is not possible ||^, because of the tendency to form a more dense 
phase. For a harmonic trap, recent computational results for T = indicate that a stable 
condensate can be formed but with no more than about 1,400 atoms for the conditions of 
the experiment j^-Q. When the number of atoms exceeds this value, it is expected that 
the cloud will collapse into a liquid or solid or into diatomic molecules. Since condensation 
of ''Li has been reported with as many as 20,000 atoms 0, there has been speculation on 
various mechanisms that might account for the observations UsHTll. 

In this report, we present a self-consistent mean field theory for trapped bosonic atoms. 
Whereas many previous theoretical studies concentrated on the condensate at T = 0, we 
report fully quantal results for T > 0, even for T > T^.. From these results, we can study the 
rapid changes in size and especially in the shape of the atomic cloud as the temperature goes 
through Tc for comparison with experimental observations. We will show that mean-field 
theory calculations with ^Li parameters are consistent with the observed spatial distribution 
of cold atoms in the trap |0. 

The Fermi-Thomas method or local density approximation (LDA) has yielded many 
useful results for BEC This approximation is an appropriate and efficient method 



for large numbers of atoms and a > 0. On the other hand, the phenomena of interest with 



a < occur on a scale equal to or less than the dimension of the lowest quantum state, so the 
applicability of the LDA approach is severely limited in this case. Other approaches to BEG 
for ^Li have involved three-body collisions [|I^ and higher order terms in the atom-atom 
interaction constructed with derivatives of 6{fi2) 

The mean field approach to BEG is based on the pseudopotential form of the atom-atom 
interaction |]l5[|l6l. We conjecture that the atom-atom interaction for densities of interest 



(an^^^ <^ 1) can be written V{fi2) = {Ajcafi^ /M) 6 {f 12) even for a < 0. From the Li2 energies 
PJl^, the next term in the scattering series, the scattering range, may be computed. We 
find that it becomes significant at mK temperatures, but our concern here will be restricted 
to T < 1/iK. 

The discussion of the present computational approach can begin with the Gross- 
Pitaevskii equation (GPE) for T = |[T7|JT8|1 . With the above form for the atom-atom 
interaction, the GPE for atoms of mass M is 



Here V{f) is the external harmonic potential and jl the chemical potential. ip is normalized 
so that / df\il)\^ = 1. For applications to the trap used for "^Li atoms, the external potential, 
V = (M/2)(u;2^2 ^ ^2-2) _ {Mujy2){Xlz^ + X^), where u = {00, + 2ujp)/3 is the average 
oscillator frequency and Aj = uji/u for i = z or p. For a convenient distance scale, we let 
f = af, with = h/{2Muj) and then let p = fl/huj. With / (ir|0(r)p = 1, and since 
df = a^dr, (j){f) = a^^'^ip^f). Thus the scaled GPE equation is 



-v' + \{xy + xy)+gNm 



0(f) = /i0(f) (2) 

The coefficient is g = Sna/a = 8na^2Muj /h. ^Li is notable not simply because a < 0, but 
also because \g\ is relatively small. For ^'^Na, \g\ is about 6 times, and for ®^Rb, \g\ is about 
14 times the value for ^Li. So in this respect, ^Li is significantly closer to an ideal gas than 
the heavier alkali atoms. 

For T > 0, one needs to obtain the wavefunctions for many orbitals self-consistently, 
analogous to the Hartree-Fock theory of atoms. Basic equations for this approach can be 



found in Ref. |19|. With the pseudopotential approximation, the Hamiltonian equation for 



each basis state in a self- consistent set is H(j)i{r) = ei(f)i{r), with 

H = -V' + \xy + Ixy + g^NjlUr-)]^ (3) 

j 

In the sum over states, the Bose distribution law applies: 

AT- = ^ (A) 

' {exp[ie,-eo)huj-f^)/kT]-l} 

where Qi is the degeneracy and the energy eigenvalue of the ith state. The chemical 
potential, fi, is determined from the condition J2i = ^• 

The distribution law, Eq. (4), is often derived from the grand canonical distribution 
based on an exchange of particles and energy between a large ensemble and a subensemble. 
This situation seems irrelevant to a gas of trapped atoms, which (after the evaporative 
cooling is turned off) ideally exchanges neither particles nor energy with the surroundings. 
We have therefore considered a microcanonical ensemble defined to include all configurations 
of harmonic oscillator quantum states consistent with a given (small range of) total energy 
and atomic number. The number of such possibilities increases very rapidly, so that with 
only 5 atoms and T = 20 nK for the trap of Ref. [2], more than 10^ configurations were found. 
The distribution obtained strongly resembled the Bose distribution with slight differences 
not likely to affect the present discussion. Since it is difficult to extend this microcanonical 
distribution analytically to thousands of atoms, we have adopted the Bose distribution. 
However, this point warrants further study. 

Solutions to (3) are obtained by two stages of iteration. For a chosen T and /i, it is 
first assumed that 0j(r) is given by harmonic oscillator eigenf unctions for each set of quan- 
tum numbers Nx{i), Ny{i), Nz{i). The mean field is computed from the sum. Since the 
Schrodinger equation for H (Eq. (3)) is nonseparable, the alternating direction implicit 



(ADI) method is used with successive scans in the x, y and z directions ^0|]. ADl iterations 
have the effect of transforming a function 4'x{x)<f)y{y)4'z{z) into a non-factorizable eigen- 
function 0(x, y, z) of H. With cylindrical coordinates, factors of p™+i/2 make numerical 



different ion near p = unreliable, so Cartesian coordinates are used. For a < 0, the atom- 
atom interactions produce a potential well that is typically smaller in spatial extent than 
the lowest orbital, so only the lowest quantum states are significantly affected by the term 
in g. Therefore, we needed to apply ADI to only the 4 to 50 lowest states. Other states are 
retained in the sum in the form of the harmonic oscillator (HO) functions. Typically states 
up to = QkT are needed, so as many as 10"^ HO states were included. 

In comparing calculations with experimental observations, it is pertinent to recall features 
of Bose distributions that differ from the more familiar Boltzmann forms. For a given T, 
there is a range of possible p values, and the spatial distribution of the ensemble will depend 
on as well as on T. This is illustrated in Fig. 1, which shows the distribution over p 
calculated for a harmonic trap with the parameters of the permanent magnet trap used for 
''Li. The frequency is is 112 Hz in the axial direction and 149 Hz in the transverse direction 
[2]. In this figure, the lighter curves are for a=0 (an ideal gas), while for the heavier curves, 
a = — 27.2ao, as for ''Li. The attractive interaction slightly narrows the distribution at large 
N. In the limit — p oo, one has effectively a Boltzmann distribution. As fi approaches 0, 
N increases and the distribution shrinks due to the onset of condensation. Thus from Fig. 

1, it is apparent that T cannot be accurately deduced from the size of an atom cloud in a 
trap unless N is also known. A similar statement applies to measurements of the width of 
the velocity distribution, since the transformed wavefunctions (f)i{p) could have been used in 
place of 0i(r) in the summation over discrete states. 

Formation of a condensate with a fixed number of atoms as T is varied is shown in Fig. 

2, for 20,000 atoms and a = -27.2 qq. The condensate peak develops over a very narrow 
range of T. As T is reduced to slightly below 162 nK, the computation becomes numerically 
unstable due to the attractive atom-atom interactions. Similarly, for ID calculations at 
T = no eigenfunction can be found beyond a certain critical value of gN. The 
maximum number of atoms in the ground state with the anisotropic trap and with the ADI 
algorithm is about 1,000, or somewhat less than obtained with an isotropic trap of the same 
mean harmonic frequency. 
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Because of this instability, the calculations indicate that for more than about 1,000 
atoms, phase-coherent BEC in a trap cannot be achieved with ^Li. However, the term 
"Bose-Einstein condensation" often refers to phenomena occurring over the entire range 



T < Tc rather than simply at T = |T^. In a harmonic trap, the signature of BEC is the 
rapid decrease of cloud size or narrowing of the velocity distribution, and these phenomena 
occur as the temperature is lowered through Tc. A general view of this phenomenon is shown 
in Fig. 3. The filled circles on this figure represent p(l/e) values taken from the distributions 
shown in the previous figure (calculated with a=-27.2 oq), while all the curves are for a=0. 
The temperature at which the transition occurs is little affected by an a of this magnitude. 
For a = 0, p(l/e) drops to a limiting value equal to the size of the ground state, but for 
a = — 27.2ao, the collapse has no limit. Determination of the ultimate minimum T for a 
given would require an investigation of fluctuation effects, which are beyond the scope of 
this study. 

The long dashed line in Fig. 3 shows that the sharp decrease in p(l/e) occurs while 
Nq/N < 10~^. Thus it does appear that with ^Li, there is a regime of or T in which a 
"condensate peak" might be observable but in which there is no threat of collapse into a 
more dense phase because the accumulation of atoms in the ground state is far short of the 
critical number. From Fig. 3, the most obvious choice for a critical temperature (Tc) is 161 
nK, which marks the limit of numerically stable solutions for a = -27.2 oq. This is close to 
the value Tc = 167 nK from semiclassical theory with the present trap parameters and 
N. 

A consideration of Fig. 1-3 might suggest that p(l/e) as a measure of cloud size ac- 
centuates the condensate. If instead one plots the average energy, (E), results for a=0 are 
shown in Fig. 3 as the short dashed line and inner left scale. (We have divided by 3k for 



easy comparison with T.) As discussed by deGroot et al. [^], when there is a finite number 
of atoms in the ground state, (E) < 3kT. In Fig. 3 there is a discontinuity in the slope 
of {E) vs. T at Tc, so for a=0, the most significant decrease of {E) /3k occurs below T^.. 
{E)/3k does not mimic the sharp variation in p(l/e) because the condensate peak involves 



relatively few atoms. 

In the experiments 0, the evolution of a fixed ensemble as temperature is lowered is not 
observed. Rather, atoms are repeatedly placed in the trap and evaporatively cooled. After 
each load and cool cycle, the total number of atoms, A^, is measured by optical absorption, 
and also the size of the cloud is measured by imaging a probe beam. The size can be stated 
in terms of p(l/e), the value of p for which the optical density is 1/e times the maximum. 
Since this value is necessarily deduced from a two-dimensional image of the cloud and thus 
includes the effects of all atoms in the line of sight, it is not exactly equivalent to the 
p(l/e) parameter extracted from calculations as above. Nevertheless, we examine in Fig. 
4 how BEC is manifest with data on and p(l/e), as defined above, for various assumed 
temperatures. These computational results are for T > and were obtained with a=0. 
Over the range of interest, the A^ vs. p(l/e) curves are quite fiat due to the condensation 
process shown in Figs. 2 and 3. To characterize the quantum state distribution, we have 
plotted a few values of Nq/N as an orthogonal grid. The data point reported in Ref. is 
clearly within the condensation regime. 

An alternative discussion of the observations on ^Li has been developed in terms of 
rotation of the cloud [Q. It has been found (and we confirm) that if all atoms have, for 
example, 3 quanta of angular momentum around the symmetry axis, about 8,000 atoms in 
the ground state will be stable. We find further that if all atoms have = 30, a condensate 
of about 50,000 atoms will be stable. However if the total value of were not an integer 
times A^, there will be differing rates of rotation, hence a possibility for collisions, in which 
case the distribution over rrii would be statistical. To consider such a distribution, we have 
added an angular velocity term to the Bose distribution law, Eq. (3), giving an exponent 
[(e^ — eQ)huj — fi — '-fm.il/kT, where mj is the value of for the ith level and 7 is an additional 
Lagrange multiplier analogous to that used for fermions in the "cranking model" of nuclear 
physics |^3[. With this distribution, one does not obtain a concentration of atoms at a single 



high value of m,. It seems possible that the cloud of atoms in the trap may randomly acquire 
angular momentum when loaded, but this is not needed to explain the observations. 



In conclusion, one may debate the semantic question as to what consitutes BEC, but self- 
consistent mean field results presented here for ''Li atoms show that a "condensate peak" in 
the spatial distribution developes rapidly for T near Tc. The experimental data P,p4| on 
vs. p(l/e) span the region of rapid decreasing p(l/e). It is predicted that over a small range 
of for a given T, a condensate peak does occur, but for slightly greater or T slightly 
less, the cloud will collapse. Confirmation of this and other predictions with experimental 
data would be useful tests of T > mean field theory for bosons. 

The author is indebted to R. Hulet, C. Bradley and C. Sackett for many communications 
regarding the experiments on ^Li, and for helpful discussions with E. Shuryak, C. N. Yang, 
N. Balazs, J. Marburger, M. Edwards, H. Metcalf and P. Koch. This work was supported 
by the NSF, ONR, and by a grant of computer time from the Cornell Theory Center. 
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FIGURES 

FIG. 1. Calculated distribution of ^Li atoms in a harmonic trap as at Rice University, with T 
= 150 nK. The plotted curves show distributions for several values of the chemical potential, 
which determines the atom number, N. The light lines are for an ideal gas with a=0, while the 
heavy lines pertain to a=-27.2 oq, as for '^Li. The curve with smallest N is identical to that for a 
Boltzmann distribution. 

FIG. 2. Calculated spatial distribution of 20,000 ^Li atoms in a trap as in Fig. 1, with 
a = — 27.2ao. The "condensate peak" arises over a small range of T. The peak includes atoms 
in the several lowest quantum states, not simply the ground state. 

FIG. 3. Condensation phenomena near T = Tc, for 20,000 ^Li atoms in a trap as in Fig. 1. The 
solid line gives the p value at 1 /e times the maximum density (whether the distribution is Gaussian 
or not) as a function of intensity for a=0 (ideal gas), while the circles give the same parameter 
calculated with a = — 27.2ao, taken from the distributions shown in Fig. 2. (Both refer to the outer 
left axis.) The long dashed line gives Nq/N on a logarithmic scale (right axis). The short dashed 
line and inner left axis show {E)/3k from averaging over all quantum states in the distribution. 

FIG. 4. Computational results obtained with a=0, plotted in the form of parameters measured 
experimentally, namely the total number of atoms N and the 1/e width of the distribution in 
the p direction. The roughly horizontal lines are for constant T as shown at the left, while the 
approximately vertical lines indicate various values of Nq/N. A data point from Ref. 2 is shown. The 
rapid variation of p{l/e) with N for a given T reflects the occurrence of Bose-Einstein condensation. 
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